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Abstract: The Linearized Auto-Localization (LAL) algorithm estimates the position of 
beacon nodes in Local Positioning Systems (LPSs), using only the distance measurements 
to a mobile node whose position is also unknown. The LAL algorithm calculates the 
inter-beacon distances, used for the estimation of the beacons' positions, from the linearized 
trilateration equations. In this paper we propose a method to estimate the propagation of the 
errors of the inter-beacon distances obtained with the LAL algorithm, based on a first order 
Taylor approximation of the equations. Since the method depends on such approximation, 
a confidence parameter r is defined to measure the reliability of the estimated error. 
Field evaluations showed that by applying this information to an improved weighted-based 
auto-localization algorithm (WLAL), the standard deviation of the inter-beacon distances 
can be improved by more than 30% on average with respect to the original LAL method. 

Keywords: auto-localization; auto-calibration; local positioning systems; differential 
sensitivity analysis; uncertainty propagation 
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1. Introduction 

There are many location-aware applications that require the estimation of the position of persons 
or objects in indoor environments. Since the Global Positioning System (GPS) is not available inside 
buildings, several localization systems have been designed to work in indoor environments, which are 
commonly known as Local Positioning Systems (LPSs) [1]. These systems require the installation of 
several nodes at fixed positions (called beacons) in the indoor environment, whose positions must be 
known in advance. The location of a mobile node (attached to the person or object to locate) can 
be obtained by the trilateration method using the measured distances between the beacons and the 
mobile node. 

The manual determination of the beacons' positions (e.g., using measuring tapes or ultrasonic/laser 
rangers) is a cumbersome error-prone method. Therefore various techniques have been proposed 
to address the problem of obtaining the position of the beacons [2-11], sometimes called the 
auto-calibration or auto-localization problem. 

Typical auto-localization solutions are based on measuring distances from the beacons to a group of 
localized nodes, or a mobile node at several known positions, and performing an inverse trilateration 
(locating the beacons using the mobile node) [2]. In [4] four different mobile node positions are used 
to locate the beacons of a 3D LPS. In [3] three nodes with known positions are required plus a group 
of nodes with unknown positions. In [5] only the relative distances between four nodes mounted on a 
mobile robot are required. These methods, however, require an external localization system to obtain the 
position of the mobile node at each static location, which is not always available in indoor environments. 

Assuming that no information on the position of any node is known, the only available data are 
the distance measurements between the beacons and the mobile node. Using this information Duff 
and Muller [6] proposed a nonlinear least-square optimization algorithm, where the objective functions 
were the distance equations and the variables the coordinates of all nodes. The initial conditions for the 
algorithm (a first position estimation of all nodes) were obtained by a trial and error method, by randomly 
generating those conditions and choosing the best solution. In [7] and [8] an Extended Kalman Filter and 
an H-infinite filter are used, respectively, to obtain the position of the beacons and the mobile node. In 
both cases an initial position estimation for the beacons is obtained using dead-reckoning information of 
the mobile node. In [9], a distance matrix is formed with the range measurements between beacons and 
the mobile node at different locations. With that information a rough approximation of the inter-beacon 
ranges is obtained using an interpolation scheme. 

In [10] a solution of the auto-localization problem, based on the linearization of the trilateration 
equations, was presented. The method, known as the Linearized Auto-Localization (LAL) algorithm, 
neither requires a trial and error approximation (i.e., randomly generated positions) nor any external 
positioning information (such as dead-reckoning data). The LAL algorithm is used when all beacons are 
located in a plane (e.g., when all the beacons are located on the ceiling) parallel to the plane containing 
the mobile node trajectory. By using these conditions, all nonlinear terms of the trilateration equations 
can be grouped in additional variables in order to obtain linear auto-localization equations. 

In the present paper we propose a method to estimate the solution's standard deviation of the LAL 
algorithm based on a first order Taylor approximation. Since the method depends on the assumption that 
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the solution has a Gaussian distribution, a confidence parameter r is defined to measure the reliability of 
the estimated standard deviation. The information obtained from the solution's precision is later used to 
further improve the LAL algorithm. 

This paper is organized as follows. Section 2 describes the LAL algorithm [10]. In Section 3 
the equations to estimate the solution's standard deviation and their reliability are proposed. These 
are evaluated by simulation in Section 4 and also experimentally on an ultrasonic 3D LPS system in 
Section 6. Section 7 presents our conclusions. 

2. The LAL Algorithm 

In [10] a new solution for the auto-localization problem of 3D LPSs was presented based on the 
distance measurements between beacons and the mobile node. While other techniques require an 
external positioning system (e.g., a second LPS or an odometric sensor) to estimate the positions of 
the mobile node, the proposed method uses only the measurements available to the LPS being used for 
mobile node location. The method is based on the linearization of the trilateration equations by grouping 
nonlinear terms in additional variables. The LAL algorithm is defined for the case where all the beacons 
are located in a plane parallel to the plane containing the mobile node trajectory. For this particular 
configuration, a solvable initial subset of three beacon nodes is obtained based on [11]. Figure 1 shows 
the auto-localization configuration for n = 3 beacon nodes N^i = {1,2, 3} and m measurement points 
Nj,j = {4, . . . , m + 3} from the mobile node path. For n = 3 beacons nodes a minimum of m = 6 
measurements points are required [11], though more measurement points can be added to improve the 
estimated solution. From now on, the measurement points of the mobile node will be referred to as 
virtual nodes. 

Figure 1. Solvable node subset composed by three beacon nodes and six virtual nodes on 



For a configuration of n = 3 beacons and m = 6 virtual nodes, the trilateration equations can 
be written as a function of two groups of distances: the unknown inter-beacon distances D u = 
{d\2, dis, d 2 3} (the variables to be estimated) and the measured distances D M = {du, c? 24 , rf 3 4, . . . , d 39 } 



a plane. 
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between beacon nodes and virtual nodes (the available data). The final equation is expressed in the linear 
form [10]: 



AX = B 



(1) 
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are composed of the distance measurements from 
the virtual nodes D M . The vector X G IR 5xl includes the unknown inter-beacon distances. Once X is 
obtained the inter-beacon distances solution D u can be calculated with: 
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For more than three beacon nodes an incremental procedure is used. The mobile node is moved 
on a plane trying to obtain at least six measurements shared by subsets of three beacon nodes. The 
goal is to obtain all possible inter-beacon solutions D u applying the linearized Equation (1) on such 
subsets. In large areas, or if the range of the nodes is limited, multiple paths can be required. Once 
enough inter-beacon distances have been obtained, a 2D localization algorithm is used to locate every 
beacon position. In the LAL method an algorithm based on the Multidimensional Scaling (MDS) is 
used to obtain all the beacons' positions [12]. The MDS algorithm requires the distance matrix which 
represents the pairwise distances between all the beacons. Since for large localization areas it is difficult 
to move the mobile node in a path where all beacons are in range, some inter-beacon distances must be 
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estimated by other means like the nearest path between those beacons. However, this estimation is a 
rough approximation of the real inter-beacons distances and it has a negative impact in the precision of 
the MDS method. To avoid this problem, the LAL algorithm uses a modification of the MDS proposed 
in [13] known as LaMSM. This method locates only subsets of fully connected beacons forming local 
maps that are later merged in a global map containing all beacons. 

3. Precision of the Auto- Localization Method 

Consider that the distance measurements D M , between beacons and virtual nodes, are given by: 

d%j d%j ~\~ £ij (V) 

where dij is the true distance between beacon i = {1, 2, 3} and virtual node j — {4, . . . , m + 3}, and 
is a zero mean Gaussian distribution with variance cr?-. An estimation of the precision of the inter-beacon 
distances D u obtained by the LAL algorithm presented in Section 2 would be very useful in the posterior 
determination of the beacons' location. However, since it is clear that the resultant error distribution of 
D u is not Gaussian, a direct evaluation of such error can be a complex task. A good estimation can 
be obtained by using a first order approximation of the nonlinear equations based on a technique called 
differential sensitivity analysis [14]. 

3.1. Differential Sensitivity Analysis 

Differential Sensitivity Analysis (DSA) is a technique used to evaluate the error of a given function 
originated by the uncertainty present on its variables. This technique uses a first order Taylor expansion 
to obtain an approximation of the function's variance. 

Consider / variables y u = g u (Z),u = {1,..,/}, each one function of p uncorrected variables 
Z = {zi, . . . , z p } normally distributed with mean z k and variance erf . The variance V(y u ) of each 
variable y u can be expressed as: 

and the covariance C(y u , y v ) between variables y u and y v : 




Equations (8) and (9) can be written in matrix form to obtain the covariance matrix C (Y) of vector 

Y = [ Vl ,..., yi] T [15]: 

c(Y) = gc{z)g T do) 

where C(Z) is the covariance matrix of Z and Q is the Jacobian matrix of G = {gi, . . . ,gi} with respect 
to vector Z. 

The approximation obtained with the DSA is only valid if the real distribution of Y is close to a 
Gaussian distribution, that is, if the higher terms of the Taylor series can be neglected. It is important to 
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verify the range of application of the method, otherwise the estimated parameters could be significantly 
different from the true ones. The range of validity of this approximation will be further discussed in the 
present section. 

In the next subsection we will apply the DSA to obtain the variance of the least squares solution of 
Equation (1). 

3.2. Error Perturbation in Least Squares Solutions 

For the estimation of the error in the least squares solution, we will use an approximation proposed by 
Stewart [16]. Similarly to the differential technique previously used, it is based on the assumption that 
the error can be modelled as a Gaussian distribution. Although in [16] this approximation was obtained 
for the case when the components of matrix A and vector B are row-wise independent, we extend it to a 
more general case where all elements of A and B can be dependent. 

Let us define the observed matrices A' = A + A and B' = B + B, where A and B represent the error 
perturbation on the true matrices A and B, originated by the noisy distance measurements. Matrices A' 
and B' are the ones used to obtain the least squares solution, since we do not know the actual values A 
and B of such matrices. Equation (1) can be rewritten as: 

(A'-A)X = B'-B (11) 

left multiplying by the pseudo-inverse A' + = (A' T A'y L A' T and rearranging the equation we obtain: 

X' = X + A' + (B - AX) (12) 

where X' = A' + B' is the calculated solution we want to evaluate. The pseudo-inverse A' + can be 
expressed as a series function of matrix A [16]: 

X' = X + A + (B - Ax) + F(B - AX) (13) 

where F represents the higher order terms of A' + . 

Finally we can simplify Equation (13) to obtain the approximation: 

X' = X + A + (B-AX) (14) 

where the step from Equation (13) to Equation (14) is based on the hypothesis that A + 3> F. Notice that 
Equation (14) is not a first order approximation of Equation (13), but rather a more restricted one. 

Using Equation (14) we can obtain the desired covariance matrix C(X') of the least squares solution 
X', which is given by: 

C{X') = C{A + {B-AX)) 

= A + C(B - AX)A +T (15) 
= A + KC(D M )K T A +T 

where the 1Z is the Jacobian matrix of R = B — AX with respect to D M . 
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3.3. Covariance Matrix of the Inter-Beacon Distances 



The covariance matrix C (D u ) of the inter-beacon distances D u is obtained by applying Equation (10) 



to the distance Equation (6): 



i\uT 



C(D U ) =HC(X')H 



(16) 



where H represents the Jacobian matrix of Equation (6) with respect to X. 

Replacing Equation (15) in Equation (16) we obtain the final variance error of the LAL method: 

c(d u ) = nA + nc(D M )n T A +T n T (17) 

Considering the simplest case, where all the measured distances D M have the same variance 
of - = a 2 M , i = {1, 2, 3} j = {4, . . . , m + 3} we have that Equation (17) can be simplified as: 
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where the effect of matrix Mdop G IR 3x3 is equivalent to the effect of the geometric dilution of precision 
(GDOP) in GPS [17]. For our case, we define the Distance Dilution of Precision (DDOP) as the 
amplification of the standard deviation a M of the measurement errors D M onto the inter-beacon distances 
D u . This amplification only depends on the relative position of the nodes (beacons and virtuals). 

For each inter-beacon distance d ik i,k — {1, 2, 3}, between beacons i and k, a respective distance 
dilution of precision ddopik can be obtained by: 
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(19) 



where is the standard deviation of the distance d^ and mnoPuv denotes the (uv) element of 
matrix M DO p- 

3.4. Reliability of the Variance Estimation 

As discussed in Section 3.1, the variance estimation of the errors of D u is only valid if its distribution 
can be approximated to a zero mean Gaussian distribution. Since such condition depends on several 
variables, a parameter to measure the validity of the estimation is necessary. For example, a second order 
Taylor series approximation could be used to evaluate higher order moments of the error distribution, 
such as the skewness and kurtosis [18]. In practice, however, we found that evaluating only the error 
perturbation on the least squares solution (Section 3.2) is enough to verify the reliability of the estimated 
variance. This is because a more restricted approximation is used for the pseudo-inverse linearization 
than for the other non-linear equations. 

The first order Taylor series of the pseudo-inverse of matrix A' = A + A can be expressed as [19]: 

A'+ ^A + - A + AA + + (A T A)- 1 A T P ± 

(20) 

^A + + F 
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where P± = I — AA + is the complementary orthogonal projection onto the orthogonal space of A. Here 
F represents the effect of noise over matrix A'. Since the approximation used in Equation (14) is valid 
when A + ^> F, a comparison between these terms is proposed in [16] by using a parameter r that meets 
the condition: 

T *irSir (21) 

where \\.\\p represents the Frobenius norm and the Stochastic norm as defined in [19]. The 
parameter r represents the ratio between the error perturbation F, associated by the noisy measurements, 
and the real value of matrix A + . As the r value increases, the perturbation became increasingly important 
and can not be disregarded. 

In Appendix 7 it is shown that, assuming the general case where all elements of A and B can be 
dependent, r can be expressed as: 



t = \Jmax{m — 1, 2 x 5}||Q2^4+| 



(22) 



where the elements q uv of matrix Q e IR are obtained from: 

m— 1 

Quv = C (aim, a>kv) (23) 
fc=i 

being C (a ku , a kv ) the covariance between elements u and v of the kth row of matrix A. The covariance 
of these elements can be calculated by applying Equation (9) to the elements of matrix A defined in 
Equation (2). 

4. Evaluation of the DSA Method 

In order to evaluate the performance of the DSA method to predict the precision of the LAL 
inter-beacon solution, a LPS was simulated based on the node configuration shown in Figure 2. The 
measurement points are distributed along a circular path under the beacons, plus one measurement point 
at its center. In the following simulations, unless otherwise stated, the parameters are: radius r = 3 m, 
number of virtual nodes m = 12 and height h = 2 m. The ranging data was generated with an 
additive white Gaussian noise with zero mean and a standard deviation of 0.01 m. All simulations were 
performed 1,000 times. 

To avoid any confusion, we here define some terms used to refer to the different types of distance 
errors present in this section. We will refer as input errors those associated with the distance 
measurements D M . These errors have a known distribution (zero mean Gaussian noise) and a 
known standard deviation. The output errors are associated with the unknown inter-beacon distances 
D u obtained using the LAL algorithm. We call simulated values those obtained by simulating the 
auto-localization problem 1,000 times in order to calculate the output error statistics. In contrast the 
estimated values are those obtained using the methods and equations defined in this paper (Equations (18) 
and (22)). Finally offline estimates are those solutions obtained with the actual values of the distance 
measurements {i.e., when we know the exact position of all nodes). The online estimates are solutions 
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obtained with noisy distance measurements and the estimated position of the nodes obtained by the 
auto-localization method. 

Figure 2. LPS configuration used for simulation composed by 3 beacons and a circular path 
with radius r, height h and m virtual nodes. 




6 5 



m+3 



The offline estimation allows us to evaluate the limits of the DSA method under ideal conditions 
{i.e., using the noise-free values of matrices A and X and the vector B). The estimation will be reliable 
as long as the output distribution of distances D u remains approximately Gaussian. As proposed in 
Section 3.4, the r value should be sufficient to verify this reliability. The online method allows us to 
estimate "on-the-fly" the standard deviation of distances D u and use this information to improve the 
beacons' localization. This method will present the same limits than the offline estimation, plus the 
effects of using the noisy matrices A' and X' and vector B' . The effect of the noisy matrices is evaluated 
by performing 1 ,000 times the online method on each test. The r value will also be sufficient to verify 
the reliability of this method. 

4.1. Standard Deviation Estimation Analysis 

To evaluate the DSA method we calculated the quadratic mean output standard deviation 
a mean = 4= a/ aj 2 + 0"i3 + o"2 3 , of the inter-beacon distances D u , obtained when shifting the center 
of the virtual nodes path of the LPS configuration (Figure 2). The resultant standard deviation maps 
obtained by simulation (Figure 3(a)) and by offline estimation (Figure 3(b)) match almost exactly. The 
advantage of the offline estimation over the simulation is that the calculation processes of the maps are 
computationally more efficient. The evaluation of the standard deviation maps can be very useful to 
analyse the ideal route path of the virtual nodes. For example, in Figure 3(b) it is shown that the ideal 
position of the center of the path is near the central point of the beacons (1.5 m,1.3 m). If the center of 
the path is far from the central point, a mean can increase up to 8 cm, that is, eight times the input standard 
deviation. 

Figure 3(c) shows the estimated standard deviation obtained by the offline and online methods 
compared to the one obtained by simulation. As expected the online estimates' errors are higher than 
the ones obtained by the offline method. Also, since the calculation of the vector X is directly related 
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with the output standard deviation (Equation (16)), the online estimates worsens as the output standard 
deviation increases. 

Figure 3. Output quadratic mean standard deviation a mean obtained when changing the center 
of the virtual nodes path: (a) obtained by simulation (1,000 times); (b) estimated offline. In 
(c) estimated a mean obtained with the offline and online methods compared to simulation. 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 
Simulated standard deviation [m] 

(C) 

4.2. Reliability of the Variance Estimation 

In Section 3.4 we proposed that the parameter r can be used to weight the validity of the estimated 
output noise. To evaluate this premise we perform a simulation changing the two parameters that affect 
the output variance: the DDOP and the input variance. To change the DDOP we use the same LPS 
configuration as before but change the radius to the values r = {0.8 m, 1.2 m, 2 m} which produce 
a ddopi2 = {7.79, 3.57, 1.52}. For the input variance a standard deviation ranging from 0 to 0.21 m 
is used. 

In Figure 4(a) the simulated and estimated output standard deviation of distance d 12 is shown with 
respect to the input standard deviation. For the online estimates the mean and the 95% and 5% percentile 
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are shown. The graph corresponds to a ddop\ 2 of 1.52. For an input noise less than 0.1 m, the 
output standard deviation increments almost linearly with the input standard deviation, as predicted by 
Equation (18). For higher values the simulation shows that the simulated standard deviation is always 
higher than the obtained from Equation (18). This happens because for high input noises the assumption 
that the error distribution of the output distances is close to a zero mean Gaussian distribution is no 
longer valid. This can be verified in Figure 4(b), where it is shown how the skewness of the distance d 12 
increases with the input noise. Since the skewness calculated for di 2 is positive, the estimated standard 
deviation will be always lower than the one obtained in the simulation. This shows that the estimated 
offline standard deviation should not be used for a high input noise. For example, a 0.1 m input noise 
limit can be chose for this particular LPS configuration. Figure 4(a) also shows, as expected, that the 
estimated online standard deviation presents a higher error with the increment of the input noise. 

Figure 4. Evaluation of d\ 2 when increasing the input standard deviation from 0 to 0.21 m. 
The used LPS configuration has a ddop\ 2 of 1.52. In (a) the simulated, estimated offline and 
online output standard deviation. For the online estimates the mean and the 95% and 5% 
percentile are shown; In (b) the skewness of the distance d 12 is shown. 




Intput standard deviation [m] Intput standard deviation [m] 



(a) (b) 

Figure 5 shows the offline and 95% percentile online standard deviation error in percentage (compared 
to the one obtained by simulation) as a function of the r value. These results are obtained for a 
ddopu = {7.79,3.57,1.52} and a standard deviation input error ranging from 0 to 0.21 m. The 
graphics are limited to values of 0 < r < 1.5 to focus within a reasonable range, since for values of 
r > 1.5 the errors are much higher than the ones shown in the graphics. 

Figure 5(a) shows the standard deviation error in percentage obtained with the offline method. This 
error is always under 10% for r values inferior to 0.7, and within this range it seems almost constant 
regardless of the value of r. For r > 0.7 the error increases with r. It is clear that, at this point, the 
effects of the skewness of the output distances begin to be appreciable. Notice that the value of r where 
the standard deviation error begins to increase seems to be independent of the ddop\ 2 value. This shows 
that the r value is sufficient to verify the reliability of Equation (18) for the offline method, regardless of 
the input errors or the geometry of the LPS system. 
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Figure 5. Calculated standard deviation errors of distance d 12 using a LPS configuration with 
ddopi2 = {7.79, 3.57, 1.52} and a standard deviation input error ranging from 0 to 0.21 m. 
The correlation between the value r and the estimated standard deviation error is shown for 
offline estimation in (a) and for online estimates (the 95% percentile) in (b). 



50 




0 0.5 1 1.5 0 0.5 1 1.5 



Estimated offline x Estimated online x 

(a) (b) 



A similar analysis used for the offline estimation can be used for the online case. Figure 5(b) shows 
the 95% percentile of the standard deviation error obtained with the online method in percentage. Since 
for every of the 1,000 iterations a slightly different value of r is obtained (due to noisy matrices) we 
grouped r in fixed intervals of 0.05 values (e.g., any value of r between 1 and 1.05 is represented in the 
graph as a value of r = 1). It can be seen that for r < 0.7 the standard deviation error is always under 
40%, though within this range the error is not constant and always increases with the value of r. This 
effect is generated by the noise perturbation over the matrices A' and B' . Still, we can use r to limit the 
maximum error obtained when estimating the output standard deviation. For example, to obtain errors 
not higher than 20% a value of t < 0.3 can be chosen. 

In this case, the value r where the standard deviation error is lower than a given limit is not totally 
independent of the dop\2 value, so it is possible that for higher values of DDOP or input noise the error 
could be higher. However, in practical applications, the values of the DDOP, the input noise, and r will 
mostly be lower than the ones here simulated. For example, the value of r obtained in the tests performed 
in Section 4.1 was always of r < 0.13 

5. Weighted Linearized Auto-Localization Algorithm (WLAL) 

In this section a new algorithm of auto-localization is presented based on the LAL algorithm and the 
error estimation techniques exposed in Section 3. As presented in [10], and reviewed in Section 2, the 
LAL algorithm calculates the distances between the unlocated beacons that are later used to obtain their 
relative position. The algorithm calculates the inter-beacon distances of every subset of three beacons, 
therefore more than one distance estimation is usually obtained between any pair of beacons. The final 
distance between two beacons is obtained by calculating the mean of all the available distances. We 
propose to rather use a weighted mean based on the online estimated standard deviations of the distances 
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using the method presented in Section 3. If for each inter-beacon distance dy we obtain k distances 
estimations, the respective weighted mean distance is obtained by: 

n 

E(dii)fc/(^y)* 

4 = —n (24) 

E 

k=l 

where is the /c-th distance cfy estimation and (<7y)fc the respective standard deviation that can be 

calculated from Equation (19). 

By using the estimated standard deviation for each possible subset, we can give more weight to 
solutions with the best subset of beacons. With the WLAL algorithm we expect an improvement of 
the resultant standard deviation of the inter-beacon distances. The expected variance &ij for each inter- 
beacon distance dij obtained by applying the weighted mean is: 



E 

k=l 

6. Evaluation on an Ultrasonic LPS 

We apply the DSA method on a real ultrasonic LPS, where the assumptions of co-planarity between 
beacons and the presence of only Gaussian noise in the measurements are not completely true. We want 
to evaluate the performance of the method under non-ideal conditions. The test is performed on the 
3DLocus system [4], shown in Figure 6(a), which is an acoustic LPS composed by n = 7 beacons 
deployed on a cell of 2.8 m x 2.8 m x 2.8 m. The calculated standard deviation of the distance 
measurements obtained with the 3DLocus is 0.23 mm, which is much more accurate than the one used 
in the simulations in Section 4. The height between the beacons and the virtual nodes is 1.4 m. 

Figure 6. 3DLocus beacons and virtual nodes configuration used in the test, (a) 3DLocus 
acoustic localization system; (b) Nodes configuration used in the test. 
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The linearized auto-localization method is compared with the inverse positioning method presented 
in [2]. Since the inverse positioning method requires the exact position of the mobile node, a Staiibli 
Unimation industrial robotic arm with a 50 /im accuracy was used for positioning such node. An 
approximately circular path with a radius of 0.8 m of m = 9 virtual nodes were used in this test, and a 
total of 100 measurements were made on each point. From these measurements we obtain 100 individual 
estimations of the positions of the beacons. The nodes' configuration used in this test is shown in 
Figure 6(b). 

6.1. Inter-Beacon Distances Estimation 

Figure 7(a) shows the estimated online and measured distance dilution of precision ddop 12 obtained 
with different subsets of beacons. The measured ddop 12 was calculated by dividing the standard deviation 
of the calculated distance d 12 by the standard deviation of the distance measurements D M . In most 
cases the estimated value is slightly lower than the actual ddopi 2 obtained in the 3DLocus, though the 
difference is always below 30%. One reason for this difference is that we are using a mean standard 
deviation of 0.23 mm for all the 3DLocus measurements while the actual value of the measurements' 
errors depends on various factors such as the distance and angle between nodes. In Figure 7(b) the 
histogram of the calculated distance c? 12 obtained using the beacons' subset {1, 2, 3} is shown. As can 
be observed, the error distribution of the calculated distance resembles a Gaussian distribution. In order 
to verify the assumption of a Gaussian error distribution on all the calculated inter-beacon distances, 
we run a Lilliefors test for normality [20]. The test established with a p-value of 0.05 that the error 
distribution on the calculated inter-beacon distances approximates a normal distribution (i.e., there is a 
5% probability that the normality is a false positive). 

Figure 7. Experimental evaluation of the estimated distance d\ 2 obtained with different 
beacons' subsets. In (a) the estimated online and measured ddop 12 \ In (b) the histogram of 
the calculated distance d\ 2 obtained using the beacons' subset {1, 2, 3}. 




(a) (b) 

In Table 1 a comparison of the output standard deviation a t obtained with the LAL algorithm and the 
standard deviation a w obtained with the WLAL algorithm are shown. In every case a w was lower than 
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o\. For example, for the distance d 12 the standard deviation obtained by the WLAL algorithm is 65% 
less than the one obtained with the LAL algorithm. Since we know that the ideal position of the virtual 
nodes path is under the center of the triangle formed by the subset of three beacons, it is clear that the 
solution obtained with the subset {1,2,3} will have a higher error than the one obtained with the subset 
{1, 2, 5}, when estimating the d\ 2 distance. On average, an improvement of 32.8% was obtained using 
the WLAL algorithm. 

Table 1. Standard deviation of the inter-beacon distances d^ obtained using the mean 07 and 
the weighted mean a w . 





cri[cm] 


cr w [cm} 




cri[cm] 


cr w [cm} 


d^ 


cTi[cm] 


cr w [cm} 


di2 


0.118 


0.041 


G?24 


0.085 


0.058 


d 3 7 


0.060 


0.047 


di3 


0.091 


0.064 


d 25 


0.094 


0.046 


G?45 


0.054 


0.041 


"14 


0.115 


0.103 


d 2 & 


0.104 


0.102 


G? 46 


0.074 


0.030 


di5 


0.090 


0.036 


dn 


0.058 


0.041 


G?47 


0.048 


0.029 


dw 


0.071 


0.049 


d 3 A 


0.116 


0.053 




0.060 


0.047 


dn 


0.053 


0.043 


d 3 b 


0.165 


0.144 


d 5 7 


0.044 


0.033 


d 23 


0.076 


0.055 


d 3 6 
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6.2. Beacons' Position Estimation 

In Figure 8 the beacons' position estimation obtained over one hundred trials with the LAL method, 
WLAL algorithm and the inverse positioning method is shown on the X-Y plane. A coordinate system 
was defined using the beacon 1 as origin and the beacon 2 as the X axis. As expected the standard 
deviation of the beacons' position is improved by using the weighted mean instead of a simple mean of 
the estimated inter-beacon distances. On average, an improvement of 22% was obtained on the standard 
deviation on axis X and Y using the WLAL algorithm. 

Figure 8 also shows bias errors between the auto-localization methods. Using the inverse positioning 
solution as the ground truth location of the beacons, a mean RMS error of the beacons' positions of 
0.76 cm and 0.92 cm is obtained with the LAL and WLAL algorithms respectively. The greatest bias 
observed is for the fifth beacon (1.3 cm). The presence of a bias between the different methods could 
be originated by several causes. First, although the cell structure of the 3DLocus is positioned parallel 
to the floor, it seems that there is a small inclination that causes a height variation on the beacons. The 
inverse positioning method shows a mean height variation of up to 8 mm between the beacons. Since the 
LAL and WLAL methods assume that the beacons are on a plane, the height difference will generate a 
discrepancy between the inverse positioning and the linearized methods. A second source of the observed 
bias could be the noisy distances used in the solution of the inter-beacon distances. In this paper we are 
only taking into account the influence of the noise using a first order approximation, however, it is known 
that the higher terms of the noise can also generate a bias [21]. 
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Figure 8. X-Y coordinates of 3DLocus beacons' positions obtained with the LAL algorithm, 
the WLAL algorithm and the inverse positioning method. The solid lines represent the 90% 
confidence ellipses of the estimated positions, (a) Beacon 2; (b) Beacon 3; (c) Beacon 4; 
(d) Beacon 5; (e) Beacon 6; (f) Beacon 7. 
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To evaluate these hypotheses, we simulated a LPS using the positions obtained with the inverse 
positioning as the real location of the beacons, but positioning them in the same plane (by setting the 
Z coordinate to be zero on all beacons) and using simulated measurements with the same standard 
deviation. The simulation showed no bias between the methods in any beacon, as can be seen for the 
fifth beacon in Figure 9, showing that the effect of the higher terms of the noise is negligible. When the 
beacons were not placed in the same plane, the simulation showed a bias between all the auto-localization 
methods. Using the beacon's altitude obtained with the inverse positioning, the simulation showed a bias 
of 5 mm for the fifth beacon. Finally, systematic range errors originated from the actual LPS can also 
add to the bias observed [22] (e.g., the range measurement is affected by the orientation between the 
beacon and the mobile node). Any error originated by this bias and also the small height variation 
between beacons can be compensated using an optimization algorithm such as the ones used in [6-8]. 
The WLAL solution could be used as a first estimation of the beacons position, since with a 1 .5 cm offset 
any optimization algorithm will easily converge to a more accurate solution. 

Figure 9. X-Y coordinates of beacon 5 obtained by simulation. The simulation does not 
presents the bias observed in the experimental test with the 3Dlocus LPS. 
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7. Conclusions 

In this paper a method to estimate the error of the solutions obtained with the Linearized 
Auto-Localization (LAL) method was presented. The method is based on the differential sensitivity 
analysis that uses a first order Taylor approximation to obtain the function's error variance. Since the 
method depends on such approximation, a confidence parameter r was defined to measure the reliability 
of the estimated error distribution. The differential sensitivity analysis showed that the standard deviation 
of the solution, obtained by the LAL method, is proportional to the standard deviation of the distance 
measurements and a matrix M dop dependent of the geometry of the localization system, which is similar 
to the geometric dilution of precision on localization problems. 

Two versions of the error estimation were evaluated: (a) the offline method that requires the location 
of all nodes and (b) the online method where all the nodes' positions are unknown. The first method is 
useful to evaluate possible path routes and strategies that maximize the accuracy of the auto-localization 
algorithm. The advantage of this method is that the calculation processes are computationally more 
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efficient than an evaluation based on a simulation process. Simulated test showed that r can be used 
to limit the maximum error obtained by the offline method regardless of the measured distances noise 
and the DDOP. The online method can be used to estimate "on-the-fly" the standard deviation of the 
inter-beacon distances obtained by the auto-localization algorithm, and use this information to improve 
its accuracy. It was also shown by simulation that the parameter r can be used to limit the maximum 
error obtained with the online method. 

Finally, a modification of the LAL method was used to evaluate the online method on an acoustic 
LPS. For the calculation of the inter-beacons distances, a weighted mean rather than a simple mean was 
proposed by using the output error statistics obtained by the online method. By using this Weighted 
Linearized Auto-Localization (WLAL) algorithm, the standard deviation of the calculated inter-beacon 
distances showed an improvement of 32.8% on average. The beacons' position estimation also presented 
a improvement of 22% of the standard deviation on axes X and Y. 
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Appendix — Calculation of the r Value 

To obtain the value of r we begin with the first order expansion of F derived from Equation (20): 

F = -A + AA + + {A T A)- 1 A T P ± (26) 
Let || F || s be the stochastic norm defined as: 



\\F\\ s =^{\\F\\l) (27) 

where E represents the expected value. The correspondent squared stochastic norm of F using 
Equation (26) is: 

||F||| =|| - A + AA + + (A T A)~ 1 A T P ± \\ 2 S (28) 
which based on [19] (Section 3.1.3) can be written as: 

||F||| =|| - A + AA + \\ 2 S + \\(A T A)- 1 A T P 1 _f s (29) 
Beginning with the first term of Equation (29) we have: 

|| - A + AA + \\l = E (\\ - A + AA + \\ 2 F 

<E(\\A + f F \\- AA + fp] (30) 



(31) 



< \\A + \\ F E (J| - AA + \\ 2 F 
The Frobenius norm of the second term in Equation (30) can be simplified as: 

E (j| - AA + \\ F ^j = E (trace [{AA + ) T AA + 
= E (trace (a +t A t AA 
= trace (a +t E (A 7 A) A 
= trace (A +T QA + ) 

where the elements q uv of matrix Q 6 IR 5x5 can be calculated by: 
q uv = E (row{A T ) u x column(A) v ^j 

= E (column(A) 7 x column(A) v ) (32) 
= E (a lu x a lv + a 2u x a 2v + • • • + a( OT -i) u x a( m _i)„) 
Since we assume that all elements of the random matrix A have mean zero, we have that: 

E (dku x a>kv) — C (dku, o>kv) (33) 
therefore Equation (32) can be simplified as: 

(m-l) 
fc=l 
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Since matrix Q is symmetric and also positive semidefinite (since it is the product of (A T A)), we can 
define Q 1 ^ 2 as the symmetric square root of Q where: 

(QVygi/2 = Q (35) 

Replacing Equation (35) in Equation (31) we obtain: 

E (|| - AA+Wl) = trace {A +T {Q 1 ' 2 )? Q 1 ' 2 A + ) 



|2 



(36) 



= IIQ 1/2 ^ + l 

So the first term of the squared stochastic norm of F Equation (29) can be written as: 

|| - A + AA + \\ 2 S <\\A + \\ 2 F \\Q l ' 2 A + \\ 2 F (37) 
A similar procedure can be used to simplify the second term of Equation (29), obtaining: 

\\(A T A)-'A T P ± \\ 2 S < \\Q 1/2 (A T Ar l \\ 2 F\\Pd 2 F (38) 

Hence: 

\\F\\ 2 S < \\A + \\ 2 F \\Q 1 / 2 A + \\ 2 F + \\Q x l 2 {A T A)- 1 \\l\\P x _\\l (39) 

Based on [16] (Section 7), from Equation (39) we obtain the expression of r as: 

II -Fll 

i ' f < y/max{m- 1,2 x 5}\\Q^A + \\ 2 = r (40) 
since A e rC^-^x 5 . From which we obtain Equation (22). 
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